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Summary 

A novel formulation termed the “integrated force method” 
(ifm) has been developed in recent years for analyzing 
structures. In this method all the internal forces are taken as 
independent variables, and the system equilibrium equations 
(ee’s) are integrated with the global compatibility conditions 
(cc’s) to form the governing set of equations. In ifm the cc’s 
are obtained from the strain formulation of St. Venant, and 
no choices of redundant load systems have to be made, in 
contrast to the standard force method (sfm). This property of 
ifm allows the generation of the governing equation to be 
automated straightforwardly, as it is in the popular stiffness 
method (sm). In this report ifm and sm are compared relative 
to the structure of their respective equations, their conditioning, 
required solution methods, overall computational require- 
ments, and convergence properties as these factors influence 
the accuracy of the results. Overall this new version of the 
force method produces more accurate results than the stiffness 
method for comparable computational cost. 


Introduction 

Solutions of structural mechanics problems must satisfy the 
appropriate equilibrium equations (ee’s) and compatibility 
conditions (cc’s) in addition to the constitutive relations (cr’s) 
describing material behavior. In what order, and to what 
extent, these three requirements are satisfied defines the 
method of analysis and the quality of the solution. This report 
concentrates on the relative roles of ee’s and cc’s in structural 
analysis. The constitutive relations, combinations of proper 
mathematical models and experimental results, here are 
presumed to be available in valid forms, even though that is 
not always the case. 

This report compares the two most fundamental approaches 
to analyzing finite element models of structures: the force 
method and the displacement method. The details of these two 
methods, and their relative characteristics, were discussed 
intensively three decades ago during the early evolution of 
computer-automated structural analysis. As is well known, the 
displacement method won out for computer automation in the 
form of the stiffness method (sm). As briefly discussed herein, 
the force method available at that time was based on the 
concept of redundant selections. It was the result of approaches 
developed for hand calculation in the precomputer era and 


proved inconvenient to automate and computationally more 
costly than the displacement method. 

A new version of the force method, described herein and 
compared with the displacement method, was introduced in 
references 1 and 2 and termed “the integrated force method” 
(ifm). It is shown in a comparison with the early force methods 
(ref. 1) that the ifm makes automation as convenient as it is 
with the displacement method and yet retains the known 
potential for superior stress-field accuracy of finite element 
models that is associated with force method solution 
techniques. Furthermore ifm provides a convenient way to 
enforce an additional constraint on a finite element model of 
a continuum, namely strain compatibility at the interelement 
boundaries. This constraint, usually not satisfied, appears to 
result in significant improvement in accuracy. 

The ifm integrates the system ee’s and the global cc’s in 
a fashion paralleling approaches in continuum mechanics (e.g. , 
the Beltrami-Michell formulation of elasticity (ref. 3)). The 
ifm is a natural way of integrating the use of ee’s and cc’s. 
In contrast, the classical force method (refs. 4 and 5) (referred 
to in this report as “the standard force method” (sfm)), 
satisfies the cc’s through the somewhat ad hoc and artificial 
concept of selected redundant internal forces. Consequently, 
ifm provides a strong motivation to reexamine the relative 
merits of the force and displacement methods within the 
context of the finite element idealization. A project was begun 
for that purpose, and its results are presented in this report. 

The primary conclusions of the comparison are as follows: 

(1) The ifm inherits from the sfm the ability to operate 
directly on stress parameters and thus to provide potentially 
more accurate stress results than does the displacement 
method. 

(2) The ifm equations for finite element discrete analysis 
form a well-conditioned system. 

(3) Discrete analysis solutions (stresses and displacements) 
obtained by ifm tend to converge to correct solutions more 
rapidly, in terms of the number of elements, than the same 
solution generated by the stiffness method. 

(4) Examples indicate that certain problems can be solved 
by ifm in less computation time than by sm. 

(5) Initial deformation problems are more elegantly treated 
by ifm than by sm. 

This research indicated that with further development the ifm 
can become a robust and versatile analysis formulation and 
a viable alternative to the popular displacement method. The 
nature of its equations makes the ifm an attractive candidate 


l 


for the inclusion of initial deformations from various sources 
(e.g., manufacturing tolerances or significant thermal effects 
and material nonlinearities). 

Historical Background 

It is of some interest to briefly revisit the historical evolution 
of the various formulations for solving structural mechanics 
problems, and specifically that of the displacement and force 
methods. 

The concepts of equilibrium of forces and compatibility of 
deformations are fundamental to analysis methods for solving 
problems in structural mechanics. There was a certain degree 
of asymmetry in the development and utilization of these two 
concepts, as described here. 

Early on, when hand calculations were used, the force 
method was favored because it resulted in a much smaller set 
of simultaneous equations, usually related to the redundant 
forces in a roof or bridge truss, than did the displacement 
method. With the appearance of computers that consideration 
lost its importance in favor of ease of automation and low 
computational cost. 

Equilibrium can be viewed as a more fundamental concept 
than compatibility. Engineers have a feel for it, perhaps 
because it was practiced consciously or subconsciously by the 
first builders of primitive human habitats in the dim past of 
human evolution, by the architects of magnificent edifices of 
biblical empires, and then by the builders of cathedrals, who 
faced the intricate equilibrium problems of Roman and Gothic 
arches and domes. If equilibrium was violated, or was 
precarious, the construction responded by tumbling down. 
Recall that the mere blast of horns is supposed to have caused 
certain walls in Jericho to collapse in a much simpler era of 
structural analysis and construction practices. 

The point is that equilibrium is such a natural concept that 
good engineers have always had a feel for it, and it was the 
guidance for most early achievements. In contrast, the more 
sophisticated and basically mathematical concept of compati- 
bility certainly was not central to the worries of these early 
builders and architects; it was not even known until math- 
ematicians defined it only a century ago. 

Rational principles to define the equilibrium conditions of 
mechanical forces had an early start with the work of 
Archimedes (287-212 B.C.) on levers and pulleys. A couple 
of millennia passed before the upsurge of rational scientific 
thought during the Renaissance brought about further 
significant theoretical developments. With the efforts of many 
scientists during the centuries that followed, the concepts of 
equilibrium and compatibility were finally developed in forms 
that eventually became useful for design calculations. Because 
it helps to understand how the current computer-automated 
analysis practices evolved, the history of this development is 
briefly reviewed. 


Renaissance geniuses like Leonardo da Vinci (1452-1519) 
and Gallileo (1564-1642) used the concept of equilibrium in 
their work, but even Gallileo, a professor of mathematics, did 
not possess the proper mathematical language to express the 
fundamental laws governing equilibrium in a continuum. That 
had to wait until the introduction of that language (i.e., 
calculus) by Newton (1642-1726) and Leibnitz (1636-1716). 
This new mathematical tool attracted many of the great 
minds of the Age of Enlightenment, who finally were in the 
possession of mathematical language to formulate correct laws 
of physics. Elasticity was introduced as a branch of 
mathematics, rich in the possible applications of the exciting 
new tool of calculus. The Bernoulli brothers and Euler were 
enthusiastic early proponents of the use of calculus in 
mechanics. Following in their footsteps many great scientists 
contributed to the early developments in elasticity. A 
fundamental contribution was made by Cauchy (1789-1857), 
who formulated the equilibrium equations both in the field and 
on the boundary for deformable bodies (refs. 3 and 6). 

The underlying principle behind the ee’s is force balance, 
which can be easily visualized. Equilibrium equations in 
general are not sufficient to solve a structural analysis problem; 
they have to be augmented by the compatibility conditions. 
In other words, ee’s are indeterminate in nature, and 
determinancy for a continuum is achieved by adding the 
compatibility conditions to them. 

The compatibility conditions in the field (which are the 
counterpart of Cauchy’s field equilibrium equations) were 
formulated in terms of strains for deformable solids by St. 
Venant (ref. 3) in 1864, decades after Cauchy’s equilibrium 
formulation. Again it took about three decades for the field 
cc’s of St. Venant to be expressed in terms of stresses by 
Beltrami and Michell in 1900 (ref. 6). 

The cc’s on the boundary, which are the counterpart of 
Cauchy’s stress boundary conditions and are henceforth referred 
to as “the boundary compatibility conditions” (bcc’s), eluded 
analysts until their recent formulation in 1986 (ref. 7). In 
contemporary elasticity, boundary indeterminancy is alleviated 
by using additional displacement boundary conditions (imposed 
on a kinematically stable structure) instead of the legitimate 
bcc’s. For discrete structures the cc’s were formulated in a 
way that would not be equivalent to the elasticity theory for 
a continuum (the strain formulation of St. Venant). 

The development of analysis methods for discrete structures 
was also accelerating during the nineteenth century for 
practical reasons. The upsurge in scientific discoveries during 
the previous centuries was closely followed by their 
exploitation for inventions and large-scale industrialization. 
The need to qualitatively predict the behavior of various 
machinery and constructions became a driving force for the 
development of analytical methods. 

One of the central problems was to analyze trusses and 
frames employed in bridges and buildings. Wooden trusses 
had been used since ancient times, and wooden bridges reached 
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spans of over 300 ft by the end of the nineteenth century. The 
first book on the analysis of bridge trusses was published in 
1847 by Whipple (ref. 8), who had impressive wooden bridges 
to his credit. When iron became available in sufficient quantity 
to construct iron bridges of various designs, accurate methods 
of calculating forces were needed in order to size structural 
members. 

The concepts of statically determinate and indeterminate 
structures were introduced. Equilibrium equations written at 
joints in terms of forces are sufficient only to calculate or 
graphically obtain member forces for statically determinate 
trusses. Two theoretical approaches (the force method and the 
displacement method) were developed for indeterminate 
structures in the second half of the nineteenth century. These 
are the foundations of the analytical methods used today. 

For an indeterminate truss there are more force unknowns 
than equations, thus the indeterminancy. Clebsch (1833-1872) 
noticed that, if ee’s are written in terms of nodal displace- 
ments, the number of equations and displacement unknowns 
is identical. With that observation the displacement method 
was born, but it was not useful because there was no practical 
way to solve the potentially large number of simultaneous 
equations by hand, except perhaps by relaxation methods, such 
as moment distribution introduced by Hardy Cross in the 
1930’s. 

A more useful method was introduced by Maxwell 
(1831-1879), who proposed cutting redundant members and 
introducing unknown redundant forces at the cuts. The 
remaining determinate structure is solved for both applied and 
redundant loads in order to obtain the internal forces and the 
relative displacements at the cuts for all the load systems. 
Because the ee’s for determinate trusses essentially represent 
a triangular system of equations, their solution is easily 
obtained even by hand calculations (refer to appendix A). Next, 
in order to reestablish discretized compatibility, the analyst 
sets up simultaneous equations that express the conditions at 
which the relative displacements due to the external loads are 
closed by the redundant loads. The solution of these equations 
yields the redundant forces, and superposition of the two sets 
of internal loads gives the final solution. This cumbersome 
method is known as the force method and referred to as “the 
standard force method” (sfm) in this report. This method 
became the analysis method of choice for generations of 
engineers because in conventional trusses and rigid frames the 
number of redundant forces, and therefore the number of 
simultaneous equations, was usually small— an overriding 
consideration before the dawn of the computer age. 

One can make the observation at this point that both the force 
and displacement methods centered around ee’s. Global 
compatibility was not dealt with at a similar conscious level 
because it is automatically satisfied in the displacement method 
and used only as an ad hoc device to augment the number of 
ee’s in the preceding formulation of the force method. 
Otherwise no parallel approach was developed that would 


satisfy the system ee’s and global cc’s simultaneously without 
recourse to the concept of redundant selection. 

General Description of Integrated 
Force Method 

All numerical solutions are approximate depending on the 
degree to which the ee’s and cc’s are satisfied or violated. 
It is a common observation that stress fields obtained by the 
heavily equilibrium-based finite element stiffness method (sm) 
in general satisfy neither the ee’s nor the cc’s in regions of 
stress concentrations or on the interelement boundaries 
(ref. 9). Of course sm provides solutions of acceptable 
accuracy to many complex problems that could hardly be 
analyzed only a couple of decades ago, the accuracy being 
dependent mostly on the choice of the finite element model. 
As will be shown, this accuracy and the efficiency of 
computation can be greatly improved by a more direct 
satisfaction of ee’s and cc’s. Recent research on cc’s led to 
the establishment of the novel formulation termed “the 
integrated force method of analysis” (refs. 1, 2, 7, and 10 
to 16). The ifm explicitly constrains the primary variables 
(which are the forces for discrete systems and the stress 
resultants for continua) to satisfy both ee’s and cc’s within 
an element and at nodes and, in addition, the cc’s on the 
interelement boundaries. The ifm thereby ensures the improved 
accuracy of the stress fields. Interelement equilibrium 
conditions, in general, are the only conditions not explicitly 
imposed by the ifm. The ifm has now been established for 
static, stability, and dynamic analyses of discrete systems and 
continua. The basic theory of ifm has been completed with 
the formulation of the variational functional for ifm (ref. 7). 
The stationary condition of the ifm variational functional yields 
all the known equations of structural mechanics along with 
the novel conditions identified as the boundary compatibility 
conditions. 

The ifm for discrete analysis, similar to sm, is independent 
of the concept of redundants and the basis determinate structure 
selection of the classical force method, referred to earlier as 
“the standard force method” (sfm) (ref. 1). The ifm for 
continuum analysis is based on ee’s and cc’s in the field and 
on the boundary. Contemporary analytical methods (continua 
or their discrete counterparts) irrespective of analysis 
methodology (sfm or sm) totally exclude the consideration of 
the boundary compatibility requirements because the bcc’s 
were not known. The ifm utilizes the bcc’s in analysis (refs. 7, 
15, and 16). Reference 1 compares the ifm to the sfm (ref. 1) 
and shows that the sfm developed for hand computation is a 
special version of the integrated force method. The sfm can 
be considered as a subset of the ifm, limited to static stress 
analysis. 

In this report the integrated force method is compared with 
the stiffness method. The methods are compared analytically 
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and numerically for finite element systems taking into 
consideration the following five criteria: 

(1) Computer automation 

(2) Solution accuracy 

(3) Stability of equation systems 

(4) Computational efforts 

(5) Versatility of methods 

As the ifm is of recent origin and only research-level software 
implementation exists, the comparison is restricted mostly to 
basic principles of the formulations. Simple numerical 
illustrations are included to clarify issues and to show relative 
performances. The basic difference between the ifm and sfm 
solution techniques is illustrated in appendix A for simple 
examples. 


The solution of equation (1) yields the n forces [F]. The m 
displacements [X] are obtained from the forces [F] by 
backsubstitution (ref. 1): 

(Xj = [J][[G] [F] + m (2a) 

where [J] is the (mxn) deformation coefficient matrix defined as 

[J] =m rows of [[S]“'] r (2b) 

Equations (1) and (2) represent the two key ifm relations for 
finite element analysis. The key equation of the stiffness 
method, expressing nodal equilibrium in terms of nodal 
displacement, has the following form: 

[K1 [X] = {?} + [P]' (3) 


Equations of Integrated Force Method 

Generation of IFM Equations 

The basic ifm equations introduced earlier (refs. 1 and 2) 
are presented here for completeness and for comparison with 
the sm equations.* A discrete or discretized structure for 
analysis can be designated as structure (n,m) where 4 ‘structure 1 ’ 
denotes type of structure (truss, frame, plate, shell, or their 
combination discretized by finite elements) and n,m are force 
and displacement degrees of freedoms (fof, do0, respectively. 
The structure (n,m) has m equilibrium equations and r = (n — m ) 
compatibility conditions. The m equilibrium equations 

[B] [F] = [P] 

and the r compatibility conditions 

[C][G] [F] = [<5R] 


are coupled to obtain the governing equations of the ifm as 


[B] 

[C][G] 


[F] 


_P 

SR 


or [S](F) = [PJ* 


( 1 ) 


where [B] is the {m X /i) equilibrium matrix, [C] is the (r x n) 
compatibility matrix, [G] is the {n x n) concatenated flexibility 
matrix, (P] is the w-component load vector, [<5R] is the r- 
component effective initial deformation vector, 


m = - [C] (j3) 0 


where [j3J 0 is the ^-component initial deformation vector, and 
[S] is the (n X n) ifm governing matrix. The matrices 
[B],[C],[G], and [S] are handed and they have full -row ranks 
of m, r, n , and ai, respectively. 


*A1! symbols are defined in appendix B. 


where [K] is the ( m X m) stiffness matrix and {P] ; represents 
equivalent loads caused by initial imperfections. 

From the ifm equations (eqs. (1) and (2)) and the stiffness 
method equation (eq. (3)) the following observations can be 
made: 

(1) In ifm the internal forces [Fj are calculated directly from 
the applied loads [P}\ In sm one has to calculate 
displacements first from loads (by using the load-displacement 
relation [K] [X] = [P], noting that loads and displacements 
have different dimensions and magnitudes) irrespective of the 
analysis requirements and then determine internal forces from 
displacements by backsubstitution, coordinate transformations, 
and differentiation or its equivalent. 

(2) The right-side vector [P]* of dimension n in ifm 
equation (1) is constructed from the w-component mechanical 
loads [Pj and the r-component effective initial deformations 
[SR]. The right-side vector of dimension m in sm includes both 
mechanical load vector [P] of dimension m and /n-component 
equivalent load vector [P] ; . The ifm load vector [PJ* is 
independent of the material characteristics and design 
parameters of the structure. The total sm load vector is a 
function of the material properties and design variables of the 
structure. The sm equivalent loads [P]* are nonzero even for 
compatible initial deformations that do not induce stresses in 
the structure. In other words the problem of initial 
deformations in the stiffness method is handled by the concept 
of equivalent loads that are nonzero even for a trivial situation 
of compatible initial deformation distribution when [SR] = [0]. 

(3) The ifm equation (eq. (1)) contains both ee’s and 
cc’s. The stiffness equation (eq. (3)) can be obtained from 
equation (1) by transforming variables and backsubstituting. 
Since equation (3) does not explicitly include the cc's, these 
equations cannot be manipulated to obtain tfm equations. 

The steps required to obtain the stiffness equation (eq. (3)) 
from the ifm equation (eq. (1)) are as follows (the derivation 
is for mechanical loads only): 

(1) Displacements are changed to deformations by using 
deformation displacement relations (ddr’s): 
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m = [B] r [Xj (4) 

(2) Deformations are then transformed to forces by using 
the deformation force relation 


Step 5: Solution of ifm equation. 

The matrix equation (1) is constructed from matrices [B], [C], 
and [G] and load [P]’, and its solution yields the forces. Dis- 
placements are obtained by backsubstituting from equations (2). 


m = [G][F] (5) 

(3) Force displacement relations are obtained from equations 
(4) and (5) as 


{FJ = [G]-'[B] r (X] (6) 

(4) The upper portion of the ee of the ifm governing equation 
(eq. (1), [B] (Fj = (Pj) is rewritten in displacements to obtain 
the stiffness equation 

[[B][G] _1 [B] 7 ] [X] = (Pj (7) 

or [K] [X] = (Pj, where [K] = [B][G] -, [B] T As mentioned 
earlier the integrated force method cannot be obtained from 
the stiffness formulation owing to the explicit absence of cc’s 
in that formulation, 

IFM Solution Procedure 

The ifm solution procedure is illustrated in appendix A for 
the example of a bar subjected to both mechanical and thermal 
loads. The prinicipal steps are as follows: 

Step 1: Assembly of the system equilibrium matrix [B]. 

The system ee matrix [B] is assembled from elemental 
equilibrium matrices by standard finite element techniques. 
The procedure to generate the elemental equilibrium matrix 
is presented in the section “Equilibrium Equations” of this 
report. 

Step 2: Generation of the global compatibility matrix [C]. 

The generation of the global compatibility matrix is 
presented in the section “Compatibility Conditions.” Since 
the deformation displacement relation (eq. (4)) utilizes the ee 
matrix [B] and the cc’s are obtained from the ddr’s, accuracies 
or errors in the equilibrium matrix are likely to be reflected 
in the compatibility matrix [C] . 


Equilibrium Equations 

The equilibrium equations, written in term of forces at the 
grid points of a finite element model, represent the vectorial 
summation of n internal forces [F] and m external loads [P]. 
The nodal ee in matrix notation gives rise to the (m x n) 
(«>m)-banded rectangular equilibrium matrix [B], which is 
independent of the material properties and design parameters 
of the indeterminate structure ( n,m ). For finite element analysis 
this matrix is assembled from elemental equilibrium matrices. 

The elemental equilibrium matrices [B] for bar and beam 
elements can be obtained from the direct force balance 
principle (ref. 4). For continuous structures, such as plates 
or shells, very few equilibrium matrices are reported in the 
literature (refs. 17 and 18). Equilibrium matrices for the plate 
flexure problem are given by Przemieniecki (ref. 17) and 
Robinson (ref. 18). Przemieniecki generates the matrix for a 
rectangular element in flexure from direct application of the 
force balance principle at the nodes. Robinson utilizes the 
concept of virtual work to derive the matrix for a rectangular- 
plate element in flexure. The procedures of Przemieniecki and 
Robinson are documented in their books (refs. 17 and 18) and 
are not repeated here. 

Energy-equivalent equilibrium matrices for finite element 
analysis can be obtained from the ifm variational functional 
(ref. 7). The procedure to be followed to generate an elemental 
equilibrium matrix from the ifm variational functional is 
illustrated next for the example of a rectangular-plate element 
in flexure. The portion of the ifm functional (ref. 7) that yields 
the equilibrium matrix [B] for a plate bending element has the 
following well-known form: 


U r 


■IK 


d 2 w 


d 2 w 


+ + H 


dx dy (8) 


where M x , M y> and are the plate bending moments and 


Step 3: Generation of the concatenated flexibility matrix [G]. 

The block diagonal flexibility matrix [G] is obtained as the 
diagonal concatenation of elemental flexibility matrices. The 
elemental flexibility matrices are obtained by discretizing 
complementary strain energy functionals according to standard 
flexibility techniques. 

Step 4: Construction of load vector [P] . 

The ifm load vector [P] * is assembled from mechanical loads 
[P] and initial deformations (j3] 0 as defined in equation (1). 


d 2 w d 2 w d 2 w 

dx 2 ’ dy 2 ’ dxdy 

represent the curvatures. The plate domain is D and the 
coordinates are (x,y). 

By appropriate choice of force and displacement functions 
the energy scalar U p can be discretized to obtain the elemental 
equilibrium matrix [B]. 

U p = (Xj r [B](Fj (9) 
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where the elemental displacement degrees of freedom are 
symbolized by [X] and the elemental force degrees of freedom 
by (Fj. 

The force fields have to satisfy two mandatory requirements: 

(1) The force fields must satisfy the homogeneous equilibrium 
equations (here the plate bending equations in the domain). 

(2) The force components F k (k = 1,2,. ..,9) must be 
independent of one another. This condition ensures the 
kinematic stability of the element. For the rectangular plate 
the force field is chosen in terms of nine independent forces: 

<F) = (F lt F 2t . . F 9 ) 


as 

M x = F { + F 2 x + F 3 y + F&y 


M y = F $ + F 6 x + F n y + F s xy f 


( 10 ) 


to* = F 9 J 

The variation of the normal moments in the field is linear, 
but the twisting moment is constant. The assumed moments 
satisfy the mandatory requirements. 

The displacement field that should satisfy the continuity 
condition (ref. 19) is chosen in terms of 12 variables to match 
the three dof (transverse nodal displacement w t and two 
rotations 0 xi and 9 yi per node i for the four nodes) and can be 
written in terms of Hermite polynomials as 


w(x f y) — Hoi(x)H 0 i(y)X { + H 0l (x)Hu (y)X 2 

+ H n (x)H 0 l (y)X 3 + H 0 l (x)H 02 (y)X 4 
+ H 0 l (x)H l 2 (y)X 5 + H n (x)H 02 (y)X 6 
^02 (-* ) Hq 2 (y ) X-j + Hq 2 (x)H\2 (>0^8 

+ 2 (** ) ^02 ( y )^9 + ^02 (jv)^io 

+ H Q 2 (x)H n (y)X ll + H n (x)H ox (y)X n 


(Ha) 

In equation (11a) the Hermite polynomials are defined as 


tfoiW = 


H m (x) = - 
H n (x) = 
H n (x) = 


x 3 - 3 a 2 x + 2 a 3 
4^~ 3 

x 3 - 3 a 2 x - 2 a 3 


4 a 3 


x 3 — ax 2 — a 2 x + a 3 


(Hb) 


4 a 1 


x 3 + ax 2 — a 2 x - a 3 


4a 2 


where X ]? X 2 X 12 are the 12 dof and a and b are the 
dimensions of the plate along the x and y directions, 
respectively. The Hermite polynomials for the y coordinate 
direction can be obtained by changing x and a to y and b, 
respectively, in equation (lib). The displacement field (eq. 
(11a)) gives rise to linear force distribution (eq. (10)) for the 
plate bending problem. 

The equilibrium matrix is obtained by substituting moments 
from equations (10) and displacements from equations (11) 
into the energy scalar given by equation (8) and integration. 
The equilibrium matrix thus obtained is presented in table I. 
The generation of the equilibrium matrix [B] from the ifm 
functional is a general procedure applicable to any other type 
of element. The matrix obtained from the functional is 
henceforth referred to as the consistent equilibrium matrix 

IB],. 

The two equilibrium matrices of Przemieniecki and 
Robinson, depicted in tables II and m, are compared with the 
ifm consistent matrix for response accuracy. In this study the 
analysis was carried out by ifm following the five steps given 
in the section “ifm Solution Procedure.” 

TABLE I. -IFM CONSISTENT EQUILIBRIUM MATRIX [B] s FOR 
RECTANGULAR-PLATE ENDING ELEMENT 


0 

b 

0 

-IF 

5 

0 

0 

a 

— 2a 2 
5 

-2 

0 

3 

0 

-F 

15 

— a 

2a 2 

5 

ab 

— 2 a 2 b 
5 

0 

b 

— ab 

-2 F 
5 

2a* 2 

5 

0 

0 

- a 2 
3 

a 3 

15 

0 

0 

b 

0 

IF 

5 

0 

0 

— a 

2a 2 

5 

2 

0 

- b 2 
3 

0 

-F 

15 

a 

— 2a 2 
5 

ab 

-2a 2 * 

5 

0 

-b 

— ab 

IF 

5 

— 2aF 
5 

0 

0 

a 2 

3 

— a 3 
15 

0 

0 

-b 

0 

-2F 

5 

0 

0 

— a 

-2a 2 

5 

-2 

0 

F- 

3 

0 

F 

15 

a 

2a 2 

5 

ab 

2 a 2 b 
5 

0 

-b 

— ab 

-2F 

5 

— lab 2 
5 

0 

0 

— a 2 
3 

-a 3 

15 

0 

0 

-b 

0 

IF 

5 

0 

0 

a 

2 a 2 
5 

2 

0 

-F 

3 

0 
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TABLE n.— PREZEMIENIECKI’S EQUILIBRIUM MATRIX 
[B] P FOR RECTANGULAR-PLATE BENDING ELEMENT 

[Dimension of element, (a,b).] 


0 

2 

0 

0 

0 

0 

0 

= 2 

-1 

b 





a 


I 

1 

0 

0 

0 
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0 

0 

0 
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0 

0 

-1 

1 

0 

0 

-2 

0 

2 

0 

0 

0 

0 

I 

b 

a 






-I 
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0 

0 

0 

0 

0 

0 

0 

0 

-1 

-I 

0 

0 

0 

0 

0 

0 

0 

0 

-2 

a 

0 

2 

b 

0 

0 

-1 

0 

0 

0 

0 

-1 

-1 

0 

0 

0 

0 

0 

1 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-2 

b 

0 

2 

a 
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0 

0 

0 

1 

-l 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

0 


IFM Consistent Matrix Versus Przemieniecki’s Matrix 

The elemental discretizations by the ifm and Przemieniecki 
approaches are identical in force and displacement degrees of 
freedoms, with fof = 9 and dof = 12 for both. This gives rise 
to identical dimensions of 12 x 9 for both matrices [B] 3 and 
[B]/>. In the ifm approach the moment variables are in units 
of moment per unit length (such as kip-inch per inch or 
kilogram-meter per meter). Przemieniecki concentrates the 
moments at the nodes, and this accounts for the dimension 
change between the elemental matrices [B^ and [B] P . The 
number of nonzero entries in the Przemieniecki matrix is 28. 
In contrast, there are 68 entries in the ifm matrix [B] s . The 
Przemieniecki matrix can be viewed as a lumped version, 
whereas the ifm matrix is consistent. The quality of the 
response was ascertained by using both the matrices [B] s and 
[B]^ to solve a finite element model of a damped square plate 
in flexure (plate material, steel; size, 40 in. (101.6 cm); and 
thickness, 0.2 in. (5.08 mm)) with a concentrated load at its 
center ( P = 1000 lb (453.59 kg)). Other analysis features in 


TABLE m.— ROBINSON’S EQUILIBRIUM MATRIX [B]„ FOR 
RET ANGULAR-PLATE BENDING ELEMENT 


[Dimension of element, (< a,b ).] 



the ifm procedure, such as the flexibility matrix and 
compatibility generation scheme, were kept identical for both 
cases. The deflection at the plate center was chosen as the 
parameter for comparison. This parameter was obtained in ifm 
by first calculating forces from equation (1) and then 
substituting [F] in equation (2). The displacement for this 
problem as given by Timoshenko (ref. 20) had the following 
value for the plate depicted in figure 1: 

= 0.4072 in. (10.343 mm) (12) 

The central displacements obtained with the two matrices 
([B] s and [B] P ) for the two finite element models are 
presented in table TV. Table TV shows that the central 
deflection obtained with the ifm matrix [B] s converged to 
Timoshenko’s series solution for the first model with four 
elements. Przemieniecki’s equilibrium matrix [B] P not only 
yielded a higher value for the central displacement for the 
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Figure I.— Clamped plate. Along clamped edges w = dw/dx = dw/dy = 0. 


TABLE IV.— IFM ANALYSIS RESULTS— RESPONSE BY 
PREZEMIENIECKI MATRIX [B] P VERSUS IFM 
CONSISTENT MATRIX [B] s 


[Timoshenko’s solution, = 0.4072 in. (10.343 mm).] 


Number of 
elements 

Przemieniecki matrix 

IFM matrix 

Central 

displacement, 

Error, 

percent 

Central 

displacement, 

Error, 

percent 

in. 

mm 

r - — 

in. 

mm 

4 

16 

0.5054 

.5509 

12.8372 

13.9929 

24.11 

35.29 

0.4083 

.4069 

10.3708 

10.3353 

0.20 

.07 


four-element model, but also did not show the tendency of 
monotonic convergence for the refined model consisting of 
16 elements. 


IFM Consistent Matrix Versus Robinson’s Matrix 

The ifm matrix bears a remarkable resemblance to 
Robinson’s matrix except for the following; 

(1) There is a change in sign between the two matrices. The 
reason is that the ifm sign convention is opposite to that of 
Robinson’s notations. 

(2) Some of the elements with coefficient 1/3 in Robinson’s 
matrix [B]* were changed to 2/5 in the ifm consistent matrix. 
This difference occurred for the 16 elements noted in table m. 
For the test example of a clamped plate, however, the 
discrepancy had a negligible effect on the solution for moments 
and displacements as shown in table V. 


When element geometry becomes complicated, it is rather 
difficult to generate the equilibrium matrix by vectorial 
summation of forces (ref. 17) or by repeated use of virtual 
influence coefficients (ref. 18). Therefore it is recommended 
that the elemental equilibrium matrices should be generated 
by following the ifm consistent approach (ref. 7). 

The equilibrium matrix \B] S corresponds to three displace- 
ment variables per node (displacement w and two rotations 
6 X and 0 y ), giving rise to 12 variables per node. This ideal- 
ization requires a displacement polynomial with 12 unknowns. 
Only 12-term Hermite polynomials are used to satisfy this 
requirement. The 12-term Hermite polynomial is known to 
cause convergency difficulties in the stiffness method. The 
effect of the 12-term Hermite polynomial on solution accuracy 
in the integrated force method was examined by using a 
displacement function given by equation (13) to obtain another 
equilibrium matrix, designated as [B] a . This displacement 
function has the following form: 

w(x,y) = a\ 4- <* 2 * + + & 4 X 2 + Qt$xy + a^y 2 

+ «7* 3 + <* 8 x 2 y + otgxy 2 + <* 1Q y 3 

+ <*i\x 3 y + «i2^> ?3 (13) 

The constants of the polynomial were related to nodal dof by 
following standard finite element techniques, and the 12 X 9 
equilibrium matrix [B] 0 was obtained from equations (8), 
(10), and (13). Few coefficients of the matrices [B] s and [B] a 
are different. 

In order to examine the effect on the solution accuracy, the 
clamped plate was analysed by using both equilibrium matrices 
[B] s and [B] a along with the appropriate flexibility and 
compatibility matrices as described in the subsection 44 ifm 
Solution Procedure.” The plate properties are defined in the 
main section “Solution Accuracy.” The central displacements 
obtained for three idealizations (model 1 has 4 elements, 
model 2 has 16, and model 3 has 36) by using the two matrices 
were as follows: 

(1) For ee matrix [B] s . The central displacements w c for 
the three idealizations were 0.2041, 0.2035, and 0.2035 in. 
(5.181, 5.169, and 5.169 mm), respectively. 


TABLE V — IFM ANALYSIS RESULTS-RESPONSE BY ROBINSON’S MATRIX [B]* VERSUS 
RESPONSE BY IFM CONSISTENT MATRIX [B], 


Number of 

Robinson’s matrix 

IFM matrix 


Ce 

displa 

ntral 

cement, 

w c 

Moment at plate 
center, 

K = My 

Central 

displacement, 

Moment at plate 
center, 

M X = My 


in. 

mm 

(kip in.)/in. 

(kg m)/m 

in. 

mm 

(kip in.)/in. 

(kg m/m) 

4 

16 

0.4083 

.4069 

10.3708 

10.3353 

0.193 

.241 

87.544 

109.317 

0.4083 

.4069 

10.3708 

10.3353 

-0.193 

-.241 

-87.544 

-109.317 
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(2) For ee matrix [B] a . The central displacements w c for 
the three idealizations were 0.2041, 0.2035, and 0.2035 in. 
(5.181, 5.169 and 5.169 mm), respectively. 

For the problem both matrices gave identical values for the 
central displacement. Moments obtained for both cases were 
also identical. 

In the stiffness method it is a known fact that the solution 
is sensitive to the choice of displacement fields (refer to eqs. 
(11) and (13)). In the ifm displacement fields do not have a 
significant effect on solution accuracy. For the plate problem 
the two different displacement fields given by equations (1 1) 
and (13) yielded identical results. This feature of ifm is further 
elaborated in the subsection “Element Level Effect.” 


Compatibility Conditions 

The compatibility conditions are constraints on strains, and 
for finite element models they are also constraints on member 
deformations (/?}. The ^-component deformation vector is 
defined as 


flexibility matrix [G] for the structure is a block diagonal 
matrix. It is obtained as the diagonal concatenation of element 
flexibility matrices. 

Generation of Compatibility Matrix [C] 

The compatibility matrix [C] is obtained by extending 
St. Venant’s theory of elasticity strain formulation to discrete 
structural mechanics (refs. 11 and 12). The procedure is 
illustrated by taking the plane stress elasticity problem as an 
example. The strain displacement relations (sdr’s) of the 
problem are 


du 

^ = — 
dx 

dv 

du dv 
Jxy ~Ty + Tx 


(1 7a) 


W = [G] (F) + (/3) 0 (14) 

where [0] is the total deformation, [G] is the flexibility matrix, 
[F] represents the member forces, and jj3) 0 is the initial 
deformation. When expressed in terms of force variables the 
cc’s involve two matrices, the flexibility matrix [G] and the 
compatibility matrix [C], and can be written as 

[C][G] (F) = [dR} (15) 


In the sdr’s three strains (e x , e y , and are expressed as 
functions of two displacements ( u and v). The compatibility 
constraint on strains is obtained by eliminating the two 
displacements from the three sdr’s, resulting in the single 
compatibility condition 


d\ d 2 Jxy __ 
dy 2 dx 2 dxdy 


(17b) 


This equation augments the ee’s in the ifm as given in 
equation (1). The flexibility matrix [G] is obtained from the 
discretization of the complementary strain energy by following 
standard techniques. For the plate element example the 
flexibility matrix [G] ? can be symbolically represented as 


U c = - [Fj 7 [G] JF) 


= | [Ml + M] - 2vM x M y + ( 1 + v)M 2 ^ dx dy 


(16) 


The two steps of St. Venant’s procedure to generate cc’s are 
as follows: 

(1) Establish the sdr’s. 

(2) Eliminate displacements from the sdr’s to obtain 
the cc’s. 

The equivalents of sdr’s in the mechanics of discrete structures 
are the deformation displacement relations (ddr’s); 
deformations [j8] of the discrete analysis are analogous to 
strains [e] of the elasticity analysis. The ddr’s can be assembled 
directly or obtained on an energy basis. 

The well-known equality relating internal strain energy and 
external work can be written for a discrete structure (n t m) in 
the following form: 


where D 1 — (Eh*/ 12), E is Young’s modulus, v is Poisson’s 
ratio, and h is the plate thickness. 

Substituting into equation (16) the moments M x , M y , and 
in terms of forces F u F 2 , F 9 as given by equation 
(10) and integrating yield the (9 x 9) flexibility matrix [G] e 
for the rectangular-plate flexure element. The generation of 
the flexibility matrix [G] is reported in the literature (refs. 4, 
17, and 18) and will not be repeated here. The system 


~ {Fj r (0) = 1 [X] r (P) (18a) 

where (X) are the nodal displacements. Equation (18a) can 
be rewritten by eliminating loads [P] in favor of forces [F] 
and using the ee ([B] [F] = [P]) to obtain the following 
relationship: 
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^Xj r [B][F) =^{F} r {/?} (18b) 

or 

“ [F] r ([B] r [X] — (/?)) = 0 (18c) 

Since the force vector [F] is not a null set, we finally obtain 
the following relation between member deformations and nodal 
displacements: 

m = [B] T {X] (19) 

The expression given by equation (19) is the general ddr 
applicable to finite element models whose ee’s can be 
symbolized as [B] [F] = [Pj. In the ddr n deformations are 
expressed in m displacements; thus there are r— (n — m) 
constraints on deformations that represent the cc’s of the 
structure ( n,m ). The cc’s are used to augment the ee’s, 
completing them to an n X n set. The r = (n — m) cc’s, which 
are obtained by eliminating m displacements from n ddr, can 
be symbolized in matrix notation as 

[C] 1(3} = [0] (20) 

In equation (20), [C] is the (r x n) compatibility matrix. It 
is a kinematic relationship and as such it is independent of 
design parameters, material properties, and external loads. The 
matrix is rectangular and banded. The maximum bandwidth of 
cc matrix [C] of an element in a finite element model depends 
on the force degrees of freedom (fof) of its neighboring 
elements. The maximum bandwidths of the compatibility 
conditions for the plate in flexure with fof = 9 shown in 
figure 2 are as follows: the maximum bandwidths (mbw) of 
compatibility conditions for the plate can be obtained as a 
function of the element location in the finite element model as 
Interior element mbw = 8 1 Zone 1 
Boundary element mbw = 54 Zone 2 
Corner element mbw = 36 Zone 3 
The generation of banded compatibility conditions is 
amenable to computer automation (refs. 11 and 12). The 
compatibility conditions of finite element analysis are 
illustrated by taking the example of a stiffened membrane 
(fig. 3). The finite element model of the structure consists of 
9 triangular membrane elements and 23 bars. Each membrane 
element has three fof; the bar has one fof. The total number 
of force variables is n = 50, consisting of 9 X 3 = 27 
membrane forces and 23 bar forces. The structure has 11 
nodes, each free node has two dof. Node 1 is fully restrained; 
node 11 is partially restrained. The number of displacement 
degrees of freedom m = (11x2-2 — 1) =19. It has 
r — (n — m) = 31 compatibility conditions. The cc’s in terms 


CBW - BANDWIDTH OF COMPATIBILITY CONDITIONS 
EBW - BANDWIDTH OF EQUILIBRIUM EQUATIONS 

FEM = FINITE ELEMENT DISPLACEMENT METHOD (12 DEGREES OF FREEDOM) 
IFM = INTEGRATED FORCE METHOD (NINE DEGREES OF FREEDOM) 



Figure 2.— Bandwidth of compatibility conditions and equilibrium equations 
for clamped plate. 

of member deformations as obtained from St. Venant’s strain 
formulation for finite element analysis of the stiffened 
membrane are as follows: 


04 + 026 = 0 

(21-1) 

01 + 024 = ^ 

(21—2) 

05 + 025 = 0 

(21-3) 

02 + 027 = 0 

(21-4) 

06 + 028 = 0 

(21-5) 

~ 026 + 029 = 0 

(21-6) 

06 + 030 = 0 

(21-7) 

010 + 031 = 0 

(21-8) 

08 + 032 = 0 

(21-9) 

07 + 033 = 0 

(21-10) 

011 4" 034 = 0 

(21-11) 

“ 032 4" 035 = 0 

(21-12) 

01! + 036 = 0 

(21-13) 
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(b) 


(a) Bar elements. 

(b) Membrane elements. 

Figure 3.— Idealized stiffened membrane structure. 


015 + 037 ~ 0 

(21-14) 

-0.10410, -O.74490 2 + 0 3 - O.831O0 5 

013 + 038 = 0 

(21-15) 

-O.71850 6 -O.13190 26 = O (21-28) 

012 + ^39 “ 0 

(21-16) 

- 0.722306 - 0.804907 + 0 9 - O.78690, o 

016 + 040 = 0 

(21-17) 

- 0.81990,, - O.12O503 2 = 0 (21-29) 

013 + 041 “ 0 

(21-18) 

-0.71860,, - O.7O920, 2 + 0,1 1040, 3 + 0 14 

017 + 045 = 0 

(21-19) 

- 0.69540,5 - O.85770, 6 = 0 (21-30) 

016 + 042 = 0 

(21-20) 

O.54660, 6 - O.1O680, 9 + O.8O180 2O 

020 + 043 = ^ 

(21-21) 

+ 0.806902, - 0.784806 + 0 47 = 0 (21-31) 

044 + 047 “ 0 

(21-22) 


021 + 046 = 0 

(21-23) 

The maximum bandwidths predicted from compatibility 

021 + 048 = 0 

(21-24) 

conditions and their actual values are shown in table VI. The 
average number of entries in cc’s is much smaller than the 

023 + 049 = 0 

(21-25) 

predictions from maximum bandwidth considerations. For the 
example of the stiffened membrane there are 31 cc’s; 27 of 

018 + 047 ~ 0 

(21-26) 

these cc’s have two entries each and the remaining 4 cc’s have 
six nonzero elements. Typically the compatibility matrix is 

022 + 050 = 0 

(21-27) 

much sparser than most other matrices of structural analysis. 


II 


TABLE VI . —BANDWIDTH OF COMPATIBILITY CONDITIONS 


Prediction from bandwidth 
consideration 

Actuals 

Maximum 

Minimum 

Maximum 

Average 

Minimum 

52 

18 

6 

a 3 

2 


a Rounded to next higher number. 


Physical Interpretation of Compatibility Conditions 

A finite element model has numerous interelement 
boundaries. From the viewpoint of continuum analysis the 
boundary compatibility conditions along these interelement 
boundaries also have to be satisfied. The conditions along the 
element interface can be symbolized as 

+ &rii = 0 (22a) 

where (&ri,6rh) represent the residue in the bcc’s of the two 
neighboring elements I and II, respectively. For the membrane 
problem the residue 5r can be represented in the following 
form (ref. 7): 


5 r = 



dx 


(N y — vN x )n x -f 


d_ 

dx 


(N x — vNy)n y 


(1 + *0 \^( N xy)n y + £( N xy )n, 


(22b) 


where N x , N y , and are the membrane stress resultants and 
n x and n y represent the direction cosines of the outward 
normal for the element interface. 

The equality constraint given by equation (22a) ensures 
interelement boundary compatibility. The correct stress fields 
should satisfy the interelement bcc’s given by equation (22a). 
Element interphase boundaries could be of complicated 
geometry; in consequence it is rather difficult to explicitly 
satisfy the interelement bbc’s given by equation (22a) apriori 
by choosing appropriate displacement or stress fields, or both. 
The satisfaction of interelement strain compatibility, and not 
just displacement continuity, is a neglected condition in the 
popular stiffness method formulation. In general, stresses 
obtained by displacement formulation along the element 
interface boundary satisfy neither equilibrium nor compatibility 
conditions. (Because of this limit ation in the stiffness-based 
finite element method, stress computation is typically avoided 
at the cardinal node points or along element interphases.) 

In the integrated force method interelement bbc’s are 
automatically enforced via the cc’s 

[C] m = [C] [G] [F] = [SR] 

Take the example of the stiffened membrane shown in figure 3. 


The cc’s between elements at interfaces are satisfied by the 
constraints given by equations (21-1) to (21-27). Take 
equation (21-12) for example. This constraint enforces 
interelement deformation continuity or bcc’s between 
membrane elements 26 and 27 along the interelement boundary 
defined by nodes 3 to 6 as depicted in figure 4(a). Such 
interelement cc’s are designated by the symbol (MM), 
membrane-to-membrane compatibility, in figure 4(a). Take 
the next cc given by equation (21-13). This constraint enforces 
deformation balance conditions between membrane 27 and 
bar 1 1 along the interface defined by nodes 5 and 6. These 
types of cc’s are symbolized by BM, bar-to-membrane 
compatibility, in figure 4(a). The total number of interphase 
cc’s (both MM and BM) is 27 as shown in figure 4(a). Besides 
these cc’s there are four bay cc’s; these enforce deformation 
balance conditions of six adjoining bars as shown in 
figure 4(b). 

In the integrated force method the equilibrium conditions 
are satisfied at the nodes [B] [F] = [P], the same as attempted 
by the stiffness method. In addition, the following cc’s are 
satisfied only in the ifm: 

(1) Membrane-to-membrane compatibility 

(2) Membrane-to-bar compatibility 

(3) Compatibility between a group of bars 

For this example the ee’s number 19 and the cc’s number 3 1 . 
The ifm satisfies all 50 equations. In contrast, the stiffness 
method is based on the 19 ee’s expressed in displacements; 
the 31 cc’s are more or less ignored. Since it is mandatory 
for the stress fields to satisfy the cc’s, their exclusion in the 
stiffness method reflects on the accuracy of the stress fields. 


Computer Automation of Formulations 

Modern structures, because of their complexity, have to be 
analyzed by digital computers. Computer automation is 
therefore an essential requirement of analysis formulation 
despite its other merits and limitations. The integrated force 
method satisfies this requirement, since all three matrices 
(equilibrium matrix [B], compatibility matrix [C], and 
flexibility matrix [G]) are amenable to automatic generation 
on a digital computer. Since the stiffness method is known 
to be amenable to computer automation, both the ifm and the 
sm satisfy this requirement. 


Solution Accuracy 

The accuracy of solutions obtained either by ifm or sm is 
of paramount importance in the choice of analysis formulation. 
During the formative period of finite element analysis Argyris 
(refs. 21 and 22) solved several problems both by the classical 
force method (sfm) and the stiffness method. Scrutiny of this 
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INTERFACE COMPATIBILITIES 

MEMBRANE TO MEMBRANE 
BAR TO MEMBRANE 



(a) 



1 3 3 

(b) 



(a) Interface. 

(b) Bay. 

Figure 4.— Compatibilities of stiffened membrane structure. 


study indicates that the sfm predicted more accurate stress 
solutions for almost all examples, including plates and 
cylinders, whereas inaccuracies can be noticed in the stiffness 
computations. The ifm retains all the favorable features of the 
classical force method (sfm) as far as solution accuracy is 
concerned. Consequently ifm predictions have to be as accurate 
as Argyris’s sfm results. 

To check the accuracy of solutions between the integrated 
force method and the stiffness method in the context of finite 
element analysis, we developed two plate bending elements 
for the ifm: 

(1) A rectangular element with four nodes 

(2) A triangular element with three nodes 

The quality of solutions obtained by the ifm and the sm was 
compared for both types of elements by taking the example 
of clamped plate bending under a concentrated load. The plate 


parameters were as follows: 

Size of plate, a = fr, in. (cm) 40 (101.6) 

Thickness of plate, h, in. (mm) 0.2 (5.08) 

Young’s modulus, E , ksi (kg/mm 2 ).... 30 000 (21 091.81) 

Poisson’s ratio, 0.3 

Magnitude of concentrated load at 
center, P, lb (kg) ..500 (226.795) 


In the stiffness method nodal stress parameters calculated by 
backsubstituting from grid point displacements are discon- 
tinuous and ambiguous (ref. 9). Calculating forces at the nodes 
is routinely avoided in the stiffness method. Because of that 
the noncontroversial nodal displacement was used in the 
comparison. It should, however, be remembered that in the 
ifm forces are the primary variables from which displacements 
are obtained by backsubstitution. 

The central deflection of the plate given by Timoshenko 
(ref. 20) is 

w c = 0.2036 in. (5.715 mm) 

For the stiffness analysis two well-known analysis codes, ask a 
(ref. 23) and msc/nastran (ref. 24), were used. The types 
of plate bending elements used were 

(1) quad-4: Rectangular element with 12 degrees of 
freedom. Both aska and msc/nastran have quad-4 elements. 

(2) trib-3 : Triangular element of the aska program with 
nine degrees of freedom. 

(3) tuba-3: Higher order triangular element of the aska 
program with 18 degrees of freedom. 

(4) tria-3: Triangular element of msc/nastran with nine 
degrees of freedom. 
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The first three elements are well known in the literature and 
are popular in practice. The quad- 4, a quadrilateral element 
with six dof per node was used here as a rectangular flat-plate 
element with three dof at grid points with appropriate 
specialization. Likewise flat-plate elements were obtained from 
the triangular elements of the general-purpose program. 

Analysis by ASKA Code 

quad-4 is a rectangular element, it has three dof per node, 
consisting of transverse displacement w and two rotations (6 X 
and dy) y and its in-plane rotation is restrained. The element 
displacement field corresponds to a cubic displacement 
function with 12 constants, trib-3 and tuba-3 of the ask a 
code are triangular elements, tuba-3 is a higher order element 
with 18 degrees of freedom. It has six degrees of freedom per 
node consisting of one displacement, two slopes, and three 
curvatures, trib-3 is a plate element with 9 degrees of 
freedom. It has three degrees of freedom per node consisting 
of one displacement w and two rotations (6 X and 6 y ). 

The results obtained for the central displacement w r by 
using aska quad-4, trib-3, and tuba-3 elements are presented 
in table VII for several finite element idealizations, shown in 
figure 5. The central displacement obtained by using the aska 
quad-4 element tended to converge to Timoshenko’s solution 
for the finite element model with 64 and 100 elements at a 
residual error of about 2 percent. As table VH shows, there was 
no substantial difference in the convergence rate for the two 
triangular elements. As before, the stiffness method required a 
very fine mesh of about 128 elements to show some convergency 
to the closed-form solution. For the 128-element model the 
residual error was about 10 percent for tuba-3 and 6 percent 
for trib-3. 

MSC/nastran Verification 

The example of a clamped plate under a central concentrated 
load as shown in figure 1 was again solved, this time by using 
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(b) 



(a) Rectangular elements. 

(b) Triangular elements. 

Figure 5.— Discretization of plate by using rectangular and triangular finite 
elements. 


the quad-4 element of the nasa structural analysis code 
msc/nastran. The quad-4 elements of the msc/nastran and 
aska codes are identical with respect to the two fundamental 
element characteristics, such as the number of nodes and the grid 
point displacement degrees of freedom. Results obtained for three 
models with 4, 16, and 36 elements, respectively, are presented 
in table VTIL 

As tables VII and Vm show, the agreement in the results 
obtained by the two codes (msc/nastran and aska) for the 
problem was good. For 4 and 16 elements the aska results 
were marginally better than the msc/nastran solution. 
However, msc/nastran showed monotonic convergence with 
a residual error of 0.29 percent for the 36-element model, 


TABLE vn.— ASKA CODE ANALYSIS RESULTS 


[Timoshenko’s solution, w c = 0.2036 in. (5.715 mm).] 


Number of 
elements 

QUAD^t rectangular element 

TRIB-3 triangular element 

TUBA-3 triangular element 

Central 

displacement, 

Error, 

percent 

Central 

displacement, 

Error, 

percent 

Central 

displacement. 

Error, 

percent 

in. 

mm 

in. 

mm 

in. 

mm 

4 

8 

16 

32 

64 

100 

128 

0.0409 

1.0388 

80.00 

0.0860 

.0650 

2.1844 

1.6510 

57.76 

68.07 

0.077! 

.0641 

1.9583 

1.628 

62.13 

68.51 

.1995 

5.0673 

2.W 

.1547 

3.9294 

24.01 

.1607 

4.0818 

21.07 

.2092 

.2092 

5.3137 

5.3137 

2.75 







2.75 

.1837 

4.6660 

9.82 

.1910 

4.8514 

6.18 
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TABLE Vffl.-MSC/NASTRAN AND IFM ANALYSIS RESULTS 


[Timoshenko’s solution, w c = 0.2036 in. (5.1715 mm).] 


Number of 
elements 

NASTRAN QUAD-4 

IFM rectangular element 

Central 

displacement, 

w c 

Error, 

percent 

Central 

displacement, 

Error, 

percent 

in. 

mm 

in. 

mm 

4 

0.027 

0.686 

86.75 

0.2041 

5.1842 

0.20 

16 

.1914 

4.862 

6.02 

.2035 

5.1689 

.05 

36 

.2042 

5.187 

.29 

.2035 

5.1689 

.05 


whereas the aska code exhibited a residual error of about 
2 percent for a 100-element model. 

Analysts by Integrated Force Method 

The ifm rectangular element assumes cubic distribution for 
displacements (eq. (11) or (13)) and linear distribution for forces 
(eq. (10)). This element corresponds to the quad-4 element, 
which also has three displacement degrees of freedom per node, 
consisting of displacement w and two slopes (6 X and 6 y ) for its 
four nodes. As far as displacement and force degrees of freedom 
and number of nodes are concerned, all three elements (ifm, 
aska, and msc/nastran) are equivalent. The results obtained 
by ifm for three different finite element models are given in table 
VIIL Remember that in ifm the forces are calculated first and 
then the desired displacements are obtained from the forces by 
backsubstitution. For the purpose of comparison, however, only 
displacement is presented, since only this variable is directly 
calculated in the stiffness method as mentioned before. 

Table Vm shows that convergence occurred for the first 
model, consisting of four elements. If symmetry is taken into 
consideration, convergence will occur for a single element. The 
residual error for the four-element model was 0.2 percent, and 
it was reduced to 0.05 percent for the second model, which had 
16 elements. Results obtained for the rectangular element by 
using ifm and the two well-known displacement-based finite 
element programs (aska and msc/nastran) are presented 
graphically in figure 6. The ifm results are barely discernible 
from Timoshenko’s solution, whereas the aska and msc/nastran 
results converge slowly. 

MacNeal and Harder (ref. 25) introduced a grading scheme 
for the evaluation of finite elements, as follows: 

A = less than 2 percent error 
B = 2 to 10 percent error 
C = 10 to 20 percent error 
D = 20 to 50 percent error 
F = greater than 50 percent error 
The comparative merits of ifm, aska, and msc/nastran 
results on the basis of residual error are given in table IX, 
which is based on MacNeaTs grading scheme. The ifm element 
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Figure 6.— Convergence analysis of clamped plate with concentrated load 
by using rectangular or quadrilateral finite elements. 

scored an “A” grade for the first model, consisting of four 
elements only, msc/nastran quad-4 required 36 elements to 
reach “A’’, and aska quad-4 was unable to achieve an “A” 
even with 100 elements. 


Results for Triangular Elements 

The convergency for the problem using triangular dements 
of the ifm and stiffness codes (aska and msc/nastran) is 
presented in figure 7. The grades secured by the elements are 
given in table IX. The ifm result was discernible from the 


TABLE IX.— MacNEAL’S GRADE CARD 
(a) Rectangular element “report card” 


Number of 
elements 

IFM 

rectangular 

MSC/NASTRAN 

QUAD-4 

ASKA 

QUAD-4 

4 

A 

F 

F 

16 

A 

B 

B 

36 

A 

A 

B 

64 

— 

A 

B 

100 

— 

— 

— 


(b) Triangular element “report card” 


Number of 
elements 

IFM 

rectangular 

MSC/NASTRAN 

TRIA-3 

ASKA 

TRIB-3 

TUBA-3 

4 

B 

F 



F 

8 

A 

D 

F 

F 

16 

— 

C 

— 

— 

32 

— 

B 

C 

D 

128 

— 

— 

B 

B 


15 
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Figure 7 —Convergence analysis of clamped plate with concentrated load 
by using triangular finite elements. 


analytical solution for the first model, with four elements, but 
even so it displayed engineering accuracy. The next model, with 
eight elements, converged to the analytical solution and also 
achieved grade “A.” 

Displacement-based triangular plate bending elements are 
known to be overly stiff and exhibit slow convergence, as can 
be seen from figure 7. Neither msc/nastran nor aska elements 
could achieve grade “A.” The msc/nastran and aska analyses 
required 32 and 128 elements, respectively, to reach a “B” 
grade. 

From solutions presented in this report and from Argyris’s 
examples it is clear that the ifm yields results that are superior 
to those obtained by the stiffness method. The reason is that the 
ifm satisfies both the ee’s and the cc’s explicitly and 
simultaneously (eq. (1» whereas the stiffness method attempts 
to satisfy equilibrium equations via displacements (eq. (3)). In 
approximate methods (such as ifm and sm) one cannot expect 


accurate response prediction when some analysis equations are 
neglected, as is the case with the stiffness method. 

Element Level Effect 

To examine the effect of different types of elements on overall 
response accuracy, we obtained the central deflection of the 
clamped plate by using the following rectangular element types. 

(1) Stiffness-based plate bending elements 

(a) aska quad- 4 element 

(b) msc/nastran quad- 4 element 

(2) iFM-based plate bending elements 

(a) Przemieniecki’s element 

(b) ifm consistent element I (displacement field defined 
by eq. (11)) 

(c) ifm consistent element II (displacement field defined 
by eq. (13)) 

The elements are equivalent with respect to displacement and 
force field assumptions, being cubic and linear, respectively, 
which gives rise to 12 degrees of freedom at the four nodes 
of the element. The elements differ with respect to assumed 
polynomials, such as defined by equations (11) and (13) for 
the ifm. The results obtained for the five types of elements 
are presented in table X. As table X shows 

(1) In the stiffness method the nature of the polynomials 
used influenced the overall response to a limited extent such 
as 86 percent versus 80 percent, 6 percent versus 2 percent, 
etc., between msc/nastran and aska elements. 

(2) In the ifm the nature of the polynomials (eqs. (11) and 
(13)) had a negligible effect on the results as can be seen for 
ifm elements I and II. 

(3) Przemieniecki’s element, which was obtained by direct 
application of the law of equilibrium, can be backcalculated 
to correspond to a rather poor displacement distribution. This 
least-sophisticated element, developed three decades ago, still 
yielded acceptable response for the difficult test problem. The 
point here is that the driver for overall solution accuracy in 
the ifm is satisfaction of the system equilibrium and global 
cc’s; the element quality or the type of interpolation poly- 
nomial used play a rather less significant role. 


TABLE X.— ELEMENT LEVEL EFFECT 


Number of 
elements 

Stiffness method 

Integrated force method 

MSC/NASTRAN 

ASKA 

Consistent 

Consistent 

Lumped [B] f 




[BL 

[BL 

(Przemieniecid) 


QUAD-4 

Rectangular element 


Error in central displacement, percent 

4 

86.75 

80.00 

0.20 

0.20 

24.11 

16 

6.02 

2.04 

.05 

.05 

35.29 

36 

,29 

— 

.05 

.05 



64 

— 

2.75 






100 

— 

2.75 

— 

— 

— 
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Conditioning of Equations of 
IFM Versus SM 

Finite element analysis requires the solution of a large 
number of simultaneous equations. Hence the stability of the 
equation is a primary criterion in the choice of analysis 
formulation. The ifm matrix [S] is not symmetric, whereas 
the sm stiffness matrix [K] is symmetric. The matrix [SI is 
of higher dimension than the stiffness matrix [K], The norms 
of the ifm upper equilibrium matrix [B] and lower cc matrix 
[G] in equation (1) differ substantially. It may therefore be 
suspected that the ifm equation that contains both ee’s and cc’s 
is an ill-conditioned system. In order to examine the issue, 
we analyzed several problems (1) by scaling the cc’s and (2) 
without any scaling at all. Scaling had no effect on the solution: 
the ifm lost only one or two precision points in 14-digit 
arithmetic for almost all problems solved. 

We further compared ifm and sm equations numerically for 
a few different types of structures. The equation stability of 
sm is governed by the parameter z, defined as 


TABLE XI . — ST ABILITY OF EQU ATION S— STIFFNESS 
METHOD VERSUS INTEGRATED FORCE METHOD 


Type of 
structure 

Stiffness method, 

Z ~ \nax^\nin 

Integrated force method, 
y ^max^min 

Truss (6,4) 

16.47 

2.96 

Truss (16,13) 

218.50 

3.17 

Truss (26,21) 

1 598.10 

7.31 

Truss (36,29) 

3 378.37 

45.89 

Frame (9,6) 

1 335.96 

3.65 

Frame (18,12) 

7 840.99 

6.96 

Frame (27,18) 

21 573.55 

10.20 

Frame (36,24) 

34 533.71 

7.48 

Plate (9,6) 

4 790.65 

18.75 

Plate (18,12) 

23 759.55 

29.89 

Plate (27,18) 

51 089.50 

56.11 

Plate (36,24) 

75 191.74 

64.56 


Sparsities of [B], [C], [S], and [K] Matrices 


^max 

^■min 


(23a) 


where X max and A min are the maximum and minimum 
eigenvalues, respectively, of the symmetric stiffness matrix 
[K]. The equation stability of ifm matrix [S] is governed by 
the parameter y defined as 


y 


c 

‘“'max 

C . 
* J min 


(23b) 


where and S^n are the maximum and minimum singular 
values, respectively, of the nonsymmetrie matrix [S]. We 
evaluated z and y for several structure types and some results 
are presented in table XI. For the examples solved, the 
eigenspace of matrix [S] was much less distorted than that of 
the stiffness matrix [K]; that is, 

e - - > > 1 or z>> y (24) 

y 

Since z is greater than y, the ifm equations were more stable 
than the sm equations. One need not be influenced by the 
simple fact that the sm stiffness matrix [K], being symmetric, 
should possess better conditioning than the ifm matrix [S] . The 
reason is that the symmetric stiffness matrix [K] is a product 
of three matrices, [K] = [B] [G] _1 [B] r , and successive 
matrix multiplications and inversions deteriorate its stability. 
The ifm matrix [S] is well conditioned. The integrated force 
method thus gave rise to a superior set of equations than did 
the stiffness method. 


Comparison of Sparsities 

Matrix sparsity and bandwidth are important parameters in 
the finite element analysis. These parameters of the compati- 
bility matrix [C] were compared with the same parameters of 
the equilibrium matrix [B] for three types of structures 
designated according to the notation introduced earlier as truss 
(101,81), frame (99,66), and plate (99,66). The structures are 
shown in figure 8. The properties of the associated symmetric 
stiffness matrix [K] were included for baseline reference. The 
parameters of the matrices are tabulated in table XII. As can 
be seen, the sparsities of the ee matrix [B] and the cc matrix 
[C] were comparable but smaller than those of the stiffness 
matrix [K] . The bandwidths of the matrices [B], [C], and [K] 
were more or less comparable. Within the bands the matrices 
[B] and [C] were sparser than the stiffness matrix [K]. The 
flexibility matrix [G] is a concatenated matrix, and its maxi- 
mum bandwidth, which depends on the number of force vari- 
ables of the elements, was much smaller than dimension n of 
matrix [S], The internal sparsity of the flexibility matrix 
enhanced the zero population of the system matrix. The 
sparsities of the matrices [S], [C], and [B] for three other 
structures are shown in table XII. 

From table XII and several other examples solved, it was 
a common observation that the governing ifm matrix [S] was 
much sparser than the equilibrium matrix [B] or the stiffness 
matrix [K]. As already noted, the compatibility matrix [C] 
appeared to be the sparsest among the structural analysis 
matrices ([B], [S], and [K]). 

Sparsity Growth During Factorization of Matrix [S] 

The ifm equation system has population along two diagonals 
as depicted in figure 9(a): ee population along the top diagonal 
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(c) 


(a) Twenty-bay truss. 

(b) Eleven-story frame. 

(c) Eleven-element cantilevered plate. 

Figure 8.— Three structures used. 


and cc population along the bottom diagonal. This 
unconventional feature of the [S] matrix was avoided by 
generating adequate numbers of ee’s followed by cc’s based 
on cc bandwidth information. The process was repeated and 
a banded [S] matrix was obtained as depicted in figure 9(b). 
Sparse matrix algebra along with a nonsymmetric matrix solver 
(la05 (ref. 26)) is used for ifm analyses. In order to avoid 
the deterioration of sparsity during solution, the symmetric 
stiffness matrix is typically factored ([K] = [L] [U]), where 
the [L] and [U] are upper and lower triangular matrices. The 
population inside the band for the symmetric stiffness matrix 
[K1 is equal to the sum of its factors ([L] 4- [U]). The question 
is, What is the sparsity increase during the factoring of the 
[S] matrix? This issue ([S] = [L][U]) was examined 
numerically for a few examples, as given in table XIII. The 
factoring was carried out by the Harwell library sparse solver 
la05. In the examples solved, the actual nonzero population 
of the sum of the triangular matrices ([L] 4 [U]) was somewhat 
less than the sparsity of the fS] matrix. The program la05, 
however, printed a minor growth in sparsity ratio. This 
deviation was not a deficiency of the structure of the [S] matrix 
but was attributed to the storage scheme used and to the degree 
of compaction in the algorithm at the time of printing. 


Comparison of Computation Time 

The ifm matrix [S] is unsymmetrical and its dimension is 
(n x n). The sm matrix [K] is symmetrical and its dimension 
is {m X m) y where n> m. From this information alone, for 
identical idealization of a structure with the same number of 
elements it can be argued that solution by ifm should be 
numerically more expensive than solution by sm. To examine 
this issue, we solved three examples (a truss, a frame, and 
a plate, each with approximately 100 degrees of freedom) for 


TABLE XII. -SPARSITY OF MATRICES [B], [C], AND [K] 


Type of 
structure 

Matrix 

[Bj 

rci 

[K] 

fB] 

[C] 

[K] 

[B] 

[C] 

[K] 

Sparsity, 

percent 


Bandwidth 



Average 

Maximum 

Truss (101,81) 

4.82 

5.94 

14.5 

8 

6 

10.91 

8 

6 

12 

Frame (99,66) 

6.85 

6.52 

17.34 

11.45 

12.45 

10.91 

12 

14 

12 

Plate (99,66) 

4.53 

4.44 

25.91 

17.18 

8.82 

17.10 

18 

16 

18 


Type of 
structure 

Sparsity ratio 

[S] 



Truss (176,141) 


0.020 

0.034 

Plate (144,27) 

.075 

.126 

.039 

Plate (36,3) 

.144 

.450 

.036 
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(a) ee and cc matrices separated. 

(b) ee and cc matrices intermingled. 

Figure 9. — Population distribution in matrix [S] . 


forces and displacements by the ifm, sm, and sfm methods 
on the same computer hardware and with the same software 
techniques. The number of elements and the finite element 
model remained the same for all three analytical methods. The 
discrete element model utilized the standard flexibility and 
stiffness matrices for bar and beam elements. For the plate 
a rectangular element with three degrees of freedom per node 
was used. These were two moments and a shear force in ifm 
and two rotations and one transverse displacement in sm. 
Computation times recorded for the three examples are given 
in table XIV. Table XIV also shows the solution time by sfm 
(refs. 1, 4, and 27 to 33) based on the "turn back lu” 
procedure of Tapsu (ref. 29) and on the lp (linear 
programming) code (ref. 12). For these examples the ifm 
required less computation time than the sm. This can be 
attributed to the following factors: 


TABLE XIII.— SPARSITY GROWTH IN LOWER AND UPPER 
TRIANGULAR FACTORS OF MATRIX [S] 


Structure 

Matrix [S] 

Triangular factors 

Sparcity ratio, 
([L] + [U])/[S], 

P/A* 

Lower, 

[L] 

Upper, 

[U] 

Number of nonzeros in matrices, 
P/A 3 

Robinson’s rectangular-plate bending element 

Plate (36,3) 
Plate (144,27) 

280/230 

2056/1553 

126/15 

1263/160 

299/152 

1290/1063 

1.045/0.661 

1.202/0.789 

IFM consistent element 

Plate (36,3) 
Plate (144,27) 

300/260 

2212/1852 

140/10 

1359/78 

166/130 

1248/1104 

1.02/0.538 

1.179/0.638 


a P - number of nonzeros as printed out by Harwell subroutine LAD5; 
A = actual number of nonzeros as counted. 


TABLE XIV.— COMPUTATION TIME— INTEGRATED 
FORCE METHOD VERSUS STIFFNESS METHOD 
AND STANDARD FORCE METHOD 


Type of 
structure 

Integrated 
force method 

Stiffness 

method 

Standard 
force method 

LU a 

LP b 

Central processing unit time, 

sec 

Truss (101,81) 

0.72 

2.5 

2.14 

3.14 

Frame (99,66) 

1.52 

2.59 

2.28 

2.09 

Plate (99,66) 

1.26 

3.1 

4.06 

3.27 


“LU: turn back LU (lower and upper triangular factors of a matrix) procedure (ref. 28). 
b LP: linear programming (ref. 11). 


(1) The sm requires a series of transformations and 
backsubstitutions (from local to global system to generate 
displacements and then from global to local system to calculate 
forces). In the ifm most of the transformations are not required. 

(2) The ifm element matrices [B] and [G] can be generated 
in closed form. 

(3) For continuum analysis, operations equivalent to 
differentiation are required to obtain stresses in the sm that 
are avoided in the ifm. 

(4) Equilibrium equations represent a very sparse system 
of equations in the ifm; but for only r = ( n—m ) forces, ee s 
essentially represent a triangular system. As such a comparison 
of m sm equations to n ifm equations is not truly accurate, 
one should compare r = (n-m) of ifm to m of sm. This is 
the traditional comparison between the classical force method 
(sfm) and sm, and it holds true for ifm in a slightly different 
scale (refer to appendix A). Thus for a truss (101 ,81) the ifm 
solution time should be more or less proportional to r — 20 
and the sm solution time to 8 1 . As the steps (1) to (4) are less 
favorable for the stiffness method, the total computation 
required was less for the ifm than for the sm. In the sfm the 
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time requirement was on the higher side because generation 
of a self-stress matrix and redundant segregation involve 
considerable computation (ref. 12). 

Buckling and Dynamic Response Analyses 

Stability and vibration response analyses are routinely 
required in industry. The standard force method cannot be 
formulated to handle such response analyses directly because 
in the sfm redundants are treated separately as loads instead 
of as part of the original structure (ref. 1). The ifm has no 
such limitations and it has been extended for eigenresponse 
analysis of both discrete systems and of continua (refs. 1,10, 
and 13). To complete the comparison and illustration of the 
ifm versus the sm, we give the stability and vibration response 
analyses of the simple three-bar truss shown in figure 10. For 
these response analyses “force mode shape” was taken as the 
primary unknown. The inertia and damping parameters were 
thus written in terms of forces as 
Inertia force: 

[M][X] = [M][J][G][F] (25) 

Damping force: 

[D](Xj = [D][J][G] [F] (26) 

where [M] and [D] are the mass and damping matrices, 
respectively. 

With the typical assumption that forces and displacements 
are harmonic in time, we can write 

(F) = (27) 


/ + i 




(a) Frequency analysis. 

(b) Stability analysis. 

Figure 10.— Response analysis of three-bar truss. 


and mass of bars were neglected. The matrices [D] ; , [S], and 
[M], for a truss are as follows: 

[0^=0 (30) 


[X] = [X] m e' w (28) 

where w is radian frequency, [F] m and (X] m are the force and 
displacement mode shapes, respectively, and i is the imaginary 
unit. 

The ifm dynamic equation can be symbolized as 


[S][F] m — io> 


[P][J][G] 

[0] 


[F]„ 


[M][J][G] 

[0] 


[F) m (29a) 


[S] = 


-1 

V2 

l 

A x 


0 

1 

-1 

a 2 


V2 

1 

]_ 

A, 


(31) 


or 


[S][FL + /w[D]y{F] m + co 2 [M]y{F] m = [0J (29b) 

Here, [D] f and [M ]j are the damping and mass matrices, 
respectively, of the ifm. 

Frequency analysis by the ifm is illustrated here for a three- 
bar truss with a mass as depicted in figure 10. Damping 


y/2L(A 3 + y/2A 2 ) L(A 3 - A { ) -y/2(A l + Vl4 2 ) ] 


[M] / = 


A e E 


-y/2LA 3 -L{A 3 +A x ) -\ (2L{A X ) 

0.0 0.0 0.0 


( 32 ) 
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where A e = (A\A 2 + AjA 3 + yl2A\A-$) and A 1, ^2, an ^ ^3 are 
the areas of members 1, 2, and 3, respectively. For = 0.68 
slug (1 kgm); A u A 2 , and A 3 = 1.0, 1.0, and 2.0 in. 2 
(645.163, 645.163, and 1290.326 mm 2 ); and E = 30 000 ksi 
(21 091.81 kg/mm 2 ) the frequencies of the structure obtained 
by solving equation (29a) were 

/, = 102.93 Hz and f 2 = 155.89 Hz (33) 

The mode shapes [F], and (F) 2 for frequencies/, and f 2 , 
respectively, were 


The stability analysis was illustrated once again by using 
the example of a three-bar truss. The stability equation for 
the three-bar truss (refs. 2 and 3) for equal bar areas of 1 in. 2 
(645.163 mm 2 ) has the following form: 



l -1 l 


r -0.996 ■) 


(39) 


(F), = < 


-0.466 


(34) 


Solution of the stability equation yielded the following results: 


1.000 J 


(1) Buckling loads 


[F} 2 = 


0.259 

0.759 

1.000 


(35) 


Displacement modes were obtained from force modes by 
backsubstituting in equation (2). The displacement modes were 


[X],= 


(X) 2 = 


1.000 

-0.318 

0.318 

1.000 


(36) 


(37) 


P\,crt = 0.89 IF, P 2 .cn = 8.248E 
(2) Force mode shapes 


[F]i = 


(F] 2 = 


1.000 

0.000 

-1.000. 

0.500 

-1.000 

0.500 


(40) 


(41) 


(42) 


Vibration analysis of the structure by using the aska code was 
also performed, and identical results were obtained. Solutions 

of eigenvalue and dynamic excitation problems are reported (3) Displacement mode shapes 
in references 10 and 13. 


Stability Analysis 

The ifm has been extended to stability analysis of structures 
(ref. 1). The stability analysis equations were obtained by 
following the usual perturbation theory . The key equation for 
stability analysis by the ifm is 

[S](F] = X[K] g [J][G](F] (38a) 


[X], = 


-1.000 

0.000 


0.000 1 

1.000 j 


(43) 


(44) 


or 

[[S] - X[S]*][F] = [0] (38b) 

where [K] g is the geometric stiffness matrix and X is the 
stability parameter. The matrix [S]^ is referred to as “the ifm 
stability matrix.” 


The stability problem was also solved with the aska code and 
identical solutions were obtained for the simple problem. 

Both the integrated force method and the stiffness method 
can handle gross response analysis; however, the classical 
force method cannot be extended to dynamic and stability 
analyses of structures. 
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Versatility of Integrated Force Method 

The integrated force method is applicable to static, dynamic, 
and stability analyses of structures idealized as 

(1) Continua (here the analysis utilizes the novel boundary 
compatibility conditions) 

(2) Discrete structures such as finite element models 
The ifm has its own variational functional from which the 
analysis equations for continua and the appropriate internal 
energy equivalent matrices of finite element analysis can be 
obtained. The ifm variational functional yields all known 
equations of structural mechanics along with the additional 
bcc’s. The ifm functional can be specialized to obtain the 
potential energy and complementary energy functionals. Like 
the stiffness method, which is known to be a versatile method, 
the ifm is also a versatile analysis formulation. 


Concluding Remarks 

The features common to the integrated force (ifm) and 
stiffness (sm) methods are as follows: 

1 . Both the ifm and the sm are amenable to computer 
automation. 

2. Both formulations can handle static, dynamic, and 
stability analyses of continua and finite element discrete 
structures. 

3. Both methods have their own variational functionals. 


In ifm all internal forces are obtained from loads in a single 
step ([S] [FJ = [Pj ). Displacements are computed from 
forces by backsubstitution. In the sm displacements are 
obtained from loads first ([K] [X] = [P]); then forces are 
computed from displacements by backsubstitution. 

The ifm bestows appropriate emphasis on equilibrium 
equations and compatibility conditions. The sm emphasizes the 
equilibrium equations. In finite element analysis the ifm 
ensures the satisfaction of element interface compatibility 
conditions. In contrast, this condition is neglected in the sm 
and is expected to be satisfied by way of mesh refinements. 

Test examples show that the ifm yields accurate responses 
for both stresses and displacements even for very coarse finite 
element models. The sm requires a much finer mesh to match 
ifm accuracy. 

The ifm equations are better conditioned than the equations 
of the stiffness formulation. The ifm governing matrix [S] is 
much sparser than the stiffness matrix [K]. The ifm requires 
fewer computations than the sm to generate an accurate 
solution. 

The integrated force method is the true force method. It is 
free from the concepts of redundants and basic determinate 
structure, which deterred automation of the classical force 
method. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, April 7, 1989 
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Appendix A 

Analysis by Integrated Force Method 


This appendix is offered here as an elementary introduction 
to the basic concepts of the integrated force method (ifm) . As 
will be shown, the ifm augments the equilibrium equations 
for indeterminate structures with constraints on element 
deformations. These constraints express the requirement that 
the structural elements fit together after deformations due to 
external loads. These constraints in the ifm are written in terms 
of member forces and represent the discretized version of the 
well-known compatibility conditions of the elasticity theory . 
This is in contrast to the standard force method (sfm), where 
the additional conditions imposed are obtained by ad hoc 
selection of redundants and cuts and then requirements are 
formulated to close the gaps at the selected cuts. 

The essential difference is that the ifm uses global 
compatibility conditions, whereas the sfm uses selected local 
conditions of deformation continuity. The sfm was developed 
in the precomputer era; therefore its goal was to generate and 
solve a small set of simultaneous equations. When it was 
interpreted for computer application, three difficulties 
emerged. One was how to automate the arbitrary selections 
of the cuts, the second was that the process was more computer 
intensive than the competing displacement method, and the 
third was the higher level of difficulty in developing a library 
of finite elements than in the case of the competing 
displacement method. The ifm overcomes all these difficulties 
for the price of an unsymmetrical set of equations, which are, 
however, computationally less intentive than the symmetrical 
equations of the displacement method. 

First a very simple, perhaps overly simple problem, a fixed 
bar shown in figure 11, is considered. For this problem the 
global compatibility conditions simply state that the sum of 
element deformations vanishes after loading. This constraint 
on member deformations can be asserted by observation for 
this problem. The compatibility condition, however, is 
developed systematically. As a more general illustration a 
simple truss example follows. 

Example 1— Fixed Bar 

The ifm analysis process is illustrated by taking the example 
of a fixed bar subjected to thermomechanical loads. The bar, 
along with its end conditions and analysis parameters, is 
depicted in figure 1 1 . The total length of the bar is 3£. It is 
idealized by three finite elements consisting of a central 
element and two boundary elements of equal length t The 
cross-sectional areas of the boundary elements A x are equal; 
the area of the central element is A 2 . The bar is made of steel; 
its Young’s modulus is E\ and its coefficient of thermal 
expansion is ot. The bar is subjected to mechanical loads P \ 
and P 2 at one- third and two-thirds of its length. The 
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(a) Fixed bar. 

(b) Discretized bar (3,2). 

Figure 11.— Analysis of bar by integrated force method. 

temperature distribution of the central element is T 2 , The 
temperatures of the boundary elements are equal to T x . 

The problem is to determine the forces and the displacements 
in the bar by the ifm. The structure is discretized by using 
three elements as shown in figure 11. The force and 
displacement degrees of freedom of the structure shown in 
figure 1 1 are as follows: 

(1) Force degrees of freedom: Each element is idealized 
by one internal force. The bar has three force degrees of 
freedom (fof = n = 3; F u F 2 , F 3 ). 

(2) Displacement degrees of freedom: The bar has two 
displacement degrees of freedom, one at each of its two free 
nodes. Its dof = m = 2; X x , X 2 . The bar for the purpose of 
analysis by the ifm is designated as “bar (3,2)” It has m - 2 
equilibrium equations and r = {n — m) = 1 compatibility 
condition. 

Equilibrium equations . — The two-system equilibrium 
equations of the bar are assembled from the three elemental 
equilibrium equations: 


Element 1 ( refer to fig. 11) 


Element 2 


[Bp = 


0 

1 



[B] (2 > = 


1 

2 



(Al) 


(A2) 
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Element 3 


is required. The concatenated flexibility matrix [G] for the 
structure has the following standard form: 


[B] (3) = 


2 

0 


-1 

1 


(A3) 


The 2x3 system equilibrium matrix [0] is obtained by 
following the standard finite element assembly technique. 


[B] 


1 -1 0 

0 1 -1 


(A4a) 


The two equilibrium equations can be written as 


1 -1 

0 1 



(A4b) 
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(A8) 


Effective initial deformations. —The right side of equation 
(A6) is zero for the case of mechanical loads. If there are 
additional prescribed deformations due to some effect (e.g., 
thermal expansion), they are accommodated in the effective 
initial deformation vector [6Rj (refer to eq. (1) of the main 
text). The effective initial deformation vector is obtained from 
initial deformations $J 0 and the compatibility matrix [C], The 
initial deformation vector has the following form: 


Compatibility conditions. — The first step in obtaining the 
compatibility conditions is to establish the deformation 
displacement relations f/3] = [B] r [X], refer to eq. (19) of the 
main text). The deformation vector [fj] of dimension 3xl 
corresponds to the three elemental expansions due to forces 
E\, F 2 , and F 3 , respectively. The deformation displacement 
relation has the following form: 
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The effective initial deformation vector [<5R] is obtained from 
the formula [SR] = - [C] f/3] 0 as 
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The two displacements are eliminated from the three ddr’s 
by simple algebra to obtain one (r = n-m = 3-2=1) 
compatibility condition. The cc in terms of deformations 
([Cl [B] = [0]) has the following explicit form: 


(5Rj = — a(27’ 1 + T 2 ) (A10) 

Take as an example the temperature distributions F, = T 0 / 2 
and T 2 = — T 0 . For this temperature profile the effective 
initial deformation vector [SR] = 0. Since [SR] is zero, such 
compatible initial deformations do not induce forces in the 
structure. This fact was known before the ifm analysis was 
begun. 

From the definition of matrices [B], [C], and [G] the final 
governing ifm equation ([S] [Fj = [P]*) is assembled as 
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The cc constrains the total elongation of the bars to zero 
(0i + 0i + 03 = 0), the same as could have been asserted 
also by observation in this case. 

From equation (A6) the compatibility matrix [C] is obtained as 

[C] = [ 1 1 1 ] (A7) 

In order to express the cc given by equation (A7) in terms 
of forces (refer to eq. (5) of the main text), the flexibility matrix 


Solution of equation (All) yields the forces from which 
displacements are calculated by backsubstitution (from eq. (2) 
of the main text). 

Numerical results. — The parameters of the example are as 
follows: 

(1) Lengths of the elements: f , = f 2 = f 3 = 10 in. 

(2) Cross-sectional areas: A, = Aj= 1 in. 2 ; A 2 = 2 in. 2 

(3) Modulus of elasticity: E = 30 000 ksi 

(4) Poisson’s ratio: v = 0.3 
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Case 1, mechanical loads only.— The pressures are 
P x = 1000 kips and P 2 = 2000 kips. The internal forces are 
F) = 1400 kips, F 2 = 400 kips, and £ 3 = - 1600 kips. The 
nodal displacements are X, = 0.4667 in. and X 2 = 0.5333 in. 

Case 2, mechanical and thermal loads.— The mechanical 
loads are as given for case 1 . The thermal loads correspond to 
the following temperature distribution: T } = 0, T 2 = 2000 °C, 
and r 3 = 0. The internal forces are F { = - 40 kips, 
F 2 = - 1040 kips, and £ 3 = - 3040 kips. The nodal 
displacements are X } = - 0.0133 in. and X 2 = 1.0133 in. 

Case 3, uniform temperature.— In this case 
T\ — T 2 = 7^ = 2000 °C. The internal forces are 
F\ = F 2 = F 3 — — 3600 kips. There is no displacement 
= X 2 = 0). 


Example 2— Five-Bar Truss 

Matrix characteristics.— The characteristics of the 
governing matrices of the ifm and the sm (ifm matrix versus 
stiffness matrix [K]) are numerically illustrated by taking the 
example of the five-bar truss shown in figure 12. The design 
parameters of the bars of the truss shown in figure 12 are as 
follows: The cross-sectional areas of the five bars are 
A x = A 3 = A 4 = A 0 and A 2 = A$ — A 0 y/2. The lengths of the 
five bars are L { - = l A - a and = L 5 = a V2. The 

Young’s moduli of the five bars are E } = E 2 = 
£3 = E 4 = £5 = E. The member forces of the five bars are 
(F u F 2y £ 3 , F 4 , £ 5 ). The member deformations of the five 
bars are (ft, jS 2 , ft, jS 4 , ft). The displacements of the two 
free nodes (fig. 12) are (X u X 2 , X 3 , X 4 ). The prescribed load 
components are (P u P 2 , £ 3 , £ 4 ). For clarity the matrix 
characteristics were examined for both determinate and 
indeterminate problems. 

Determinate problem.— A four-bar determinate structure 
was generated from the indeterminate truss by eliminating the 
fifth bar. For determinate structures the equilibrium equations 
represent both necessary and sufficient conditions for stress 
analysis. There are no compatibility conditions for determinate 
structures. The ifm governing equations become 

[S] (F) = [B] [F] = [P] 


The ifm equations are obtained from the nodal force balance 
condition between the member forces [F] and the external loads 
[P] (refer to fig. 12). 
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In the displacement method the system equations are assembled 
from the elemental equations by following standard finite 
element assembly techniques. The assembled stiffness equation 
([K] [X] = [P]) for the problem has the following form: 


p 2 ,x 2 




(b) 


(a) Indeterminate truss (5,4). 

(b) Determinate truss (4,4). 

Figure 12.— Analysis of truss by integrated force method. 
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From the ifm matrix [S] (eq. (A 12)) and the stiffness matrix 
[K] (eq. (A 13)) we made the following observations: 

(1) Both are square matrices of dimension 4x4. 

(2) The matrix [S] for the determinate truss is a triangular 
matrix. For general determinate structures irrespective of their 
complexity, the governing matrix [S] after some rearrangement 
of rows and columns can be represented as a triangular matrix. 
The stiffness matrix [K] is not a triangular matrix for determin- 
ate systems as can be seen even for this trivial example. The 
triangular system of equations requires insignificant computa- 
tion, and it can be solved even manually irrespective of the 
size or complexity of the problem. This feature of the force 
method made it the popular analytical method in the 
precomputer era. The classical force method utilizes this 
property of the governing matrix even for indeterminate 
structures and reduces the number of simultaneous equations 
to the order of only the redundants. 

(3) The coefficients of matrix [S] tend to be of the same 
order of magnitude, and they are dimensionless numbers. 
These two features of the matrix ensure that matrix [S] is 
numerically stable. 

The coefficients of the stiffness matrix have the dimensions 
of force per unit length. Since the coefficients depend on 
material properties and design parameters, an ill-conditioned 
stiffness matrix can be easily obtained by changing those 
properties. 

For this simple example or any other complex determinate 
system the force method should be followed in analyses for 
the following reasons: 

(1) The most important variables, namely forces, required 
by design engineers are obtained directly in the force method. 

(2) The force methodTequfrenittTe computation; unlike for 
the displacement method there are no simultaneous equations 
to solve. 

(3) Its equations are well conditioned. 

From forces (FJ the required number of displacements can be 
obtained in a small fraction of the total analysis time. It is 
inefficient to use the stiffness method or any other formulation 
for the stress analysis of determinate problems. 


Compatibility conditions for determinate structures 
Determinate structures do not have any compatibility 
conditions (cc’s). This fact is numerically illustrated by the 
example of the determinate truss. Generation of cc’s for a 
determinate system can be attempted by following the two- 
step procedure given in the text in the section “Generation 
of Compatibility Matrix [C].” In the first step the ddr’s are 
established. In the second step the cc’s are generated from 
the ddr’s by eliminating the displacements. 

The ddr’s of the determinate structure can be obtained by 
using the equilibrium matrix [ (3 ] as (/?] - [B] r fX]. The ddr’s 
for the determinate truss have the following form: 
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(A 14) 


In equation (A 14) the bar deformations (elongations or 
contractions) are represented by (ft, ft, ft, and ft). The four 
deformations can be uniquely determined in terms of the four 
displacements as a solution to equation (A 14). In other words 
there is no constraint on deformations, and as a result deter- 
minate structures do not have any compatibility conditions. 

Indeterminate problem. — The ifm equilibrium equations 
(ee’s) for the five-bar indeterminate truss can be written as 
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(A 15) 


The ee’s of the indeterminate problem given by equation (A 15) 
still contain a triangular factor (as was the case for determinate 
structures). However, there are fewer ee’s than unknowns, 
since there are four simultaneous equations and five unknown 
forces. Equation (A 15) cannot be solved for the five unknown 
forces, hence the indeterminancy. An additional equation is 
required to augment the system to five equations in five 
unknowns, which can be solved for the five member forces. 
This additional equation is the compatibility condition. 

The cc to augment the 4 X 5 ee for the truss is obtained 
in two steps. In the first step the cc in deformations is generated 
from the ddr. In the next step the cc in deformations is written 
in forces by using the constitutive relations. 
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The ddr for the five-bar truss has the following form: 
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In equation (A 16) five deformations are expressed in terms 
of four displacements: thus there is a single (n — m = 
5 — 4 = 1) constraint on deformations that represent the cc. 
The cc is obtained by eliminating the four displacements (. X { , 
X 2 , X 3 , and X 4 ) from the five ddr’s given by equation (A 16). 
The cc has the following explicit form: 
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The cc (eq. (A 17)) of deformations is expressed in terms 
of forces by using the force deformation relations of the truss. 
These relations ([/?] = [G] (Fj) can be written in matrix notation 
as 
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Equations (A 17) and (A 18) normalized with respect to the 
flexibility parameter Li AE yield the cc of the integrated force 
method expressed in terms of force variables as 
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The cc given by equation (A 19) is added to the ee given ^ 
equation (A 15) to obtain the governing ifm equation 
«S] [F] = [P] *) as 
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Solution of equation (A20) of dimension 5x5 yields the 
five member forces F u F 2 , F 3 , F 4 , and F 5 . 

The governing equation of the stiffness method ([K] 
[X] = (Pj) for the five-bar truss has the following form: 
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The solution of equation (A21) yields the nodal displacements. 

The three conclusions for determinate structures stated 
earlier hold true even for indeterminate structures with the 
following exceptions: 

(1) The ifm yields n x n (here 5x5) equations where the 
dimension of the stiffness equation is m x m (here 4x4). 

(2) The triangular factor of matrix [S] is disturbed because 
of the compatibility conditions. 

It is common observation that 

(1) In the stiffness method compatibility is presumed to be 
satisfied apriori by appropriate choice of displacement fields, 
and the equilibrium equations ([K] [X] = [P]) are the governing 
equations. 

(2) In the classical force method (Airy’s stress function 
formulation being the typical example) the equilibrium is 
satisfied apriori (by suitable assumptions on stresses) and the 
cc’s are the governing equations. 


27 



The question is, among ee and cc which one is assumed 
and which one is enforced in the integrated force method? In 
the analysis of indeterminate structures (refer to ifm governing 
equation [S] (F) = [P] * (eq. (A20) for the indeterminate truss) 
both the ee’s and the cc’s are simultaneously and explicitly 
enforced by the integrated force method. It can be verified 
that among the five available methods of structural analysis, 
(1) the classical Airy force method, (2) the stiffness method, 
(3) the mixed method, (4) the total formulation, and (5) the 
integrated force method, only the integrated force method 
enforces both ee’s and cc’s simultaneously. 


The primary ifm variables are the forces, not a combination 
of forces, displacements, and deformations. Therefore the 
integrated force method is a direct formulation, and it should 
not be confused with the mixed method (in which forces and 
displacements are considered as m + n simultaneous primary 
unknowns) or with the total formulation (in which forces, 
deformations, and displacements are treated as the m + 2n 
primary unknowns). 
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Appendix B 
Symbols 


Ai,A 2 ,Aj 

cross-sectional areas of three truss bars 

[p] f 

equivalent loads due to initial imperfections 

a,b 

dimension of rectangular element 

(pf 

ifm load vector of dimension n 

[B] 

equilibrium matrix 

r 

number of compatibility conditions. 

[B ] P 

Prezemieniecki’s equilibrium matrix 


r = n — m 

[B] r 

Robinson’s equilibrium matrix 

[S] 

ifm governing matrix of dimension ( n X n) 

[B]j,[B] a 

ifm equilibrium matrices 

[S] b 

ifm stability matrix 

[C] 

compatibility matrix 

t 

time parameter 

[D] 

damping matrix 

Uc 

complementary strain energy 

[D]/ 

ifm damping matrix 

Up 

strain energy 

E 

Young’s modulus 

u,v 

membrane displacement components 

e 

eigen ratio 

w 

transverse displacement 

F 

internal force 

W c 

transverse displacement at center 

IF] 

force vector of dimension n 

X x ,X 2 ,...,X n 

displacement variables 

[F) m 

force mode shape of dimension n 

(X] 

displacement vector of dimension m 

/ 

frequency 

(XL 

displacement mode shape of dimension m 

[G] 

flexibility matrix of dimension ( n X n) 

x,y,z 

Cartesian coordinates 

[G) e 

element flexibility matrix 

a 

coefficient of thermal expansion 


Hermite polynomials (ij = 0, 1,2) 

m 

deformation vector of dimension n 

h 

plate thickness 

Pk 

K 111 deformation component 

i 

imaginary unity 

Mo 

initial deformation vector 

[J] 

deformation coefficient matrix of dimension 

y 

shear strain 


(m X m) 

m 

effective initial deformation vector 

[K] 

stiffness matrix of dimension ( m X m) 

Pr 

residue in bcc 

[K], 

geometric stiffness matrix of dimension 

y X y) 

strain components 


(m X n) 

6 x ,0y 

rotation about x and y axes 

L,l 

length parameters 

X 

eigenvalue 

U...U 

lengths of five bars 

V 

Poisson’s ratio 

M x ,My,M X y 

plate bending moments 

CO 

radian frequency 

[M] 

mass matrix of dimension (m X m) 



[M]/ 

ifm mass matrix 

Subscripts: 


m 

displacement degrees of freedom 



m o 

lumped mass 

crt 

critical load 

N x ,N y ,N xy 

membrane stress resultants 

max, min 

maximum or minimum values 

n 

force degrees of freedom 



n x ,n y 

direction cosines 

Superscript: 


IP] 

mechanical load vector of dimension m 

T 

transpose of matrix or vector 
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